% This is the code for calculating the magnetization dynamics of MnF2
% under a magnetic field along z axis, and a microwave is applied in xy
% plane. Here the phase of the microwave chnages from 0 to 2pi.
%parameters
pi=3.1415926;
h=6.626e-34/(2*pi);
e=-1.6e-19;
gamae=2.21*1e5/(4*pi*1e-7);
gama=gamae;
Ms=4.47E4; 
Ku=2.03E4;
Hex=53;  %exchange field
Hk=0.85;  %anisotropy field
dt=5e-15; % time interval
miu0=1;
at=10;
n=50000;
g11=1; %parameters for calculaitng spin current
g22=1;
g12=1;
g21=1;
   fh=240; %frequency
   w=fh*2*pi*1e9;   
for hn=97:1:97 
    H0=hn*0.1  % strength of the magnetic field
    alpha=0.0005; %damping constant
%H0, magnetic field
thetaH=3*3.14/180; % tilting angle from z axis
H0x=0;
H0y=H0*sin(thetaH);
H0z=H0*cos(thetaH);  

for thetaxn=1:1:37 % the phase of the microwave
%m0 intial magnetization
m1(1,1)=mfinal(hn,1);
m1(1,2)=mfinal(hn,2);
m1(1,3)=mfinal(hn,3);
m1(1,4)=1;
m2(1,1)=mfinal(hn,4);
m2(1,2)=mfinal(hn,5);
m2(1,3)=mfinal(hn,6);
m2(1,4)=1;
%Hp, setting the microwave
Hp0=0.05;
thetax=(thetaxn-1)*10*3.14/180;
thetay=0;
Hp(1,1)=Hp0*cos(thetax);
Hp(1,2)=Hp0*sin(thetay);
Hp(1,3)=0;
%exchange field
Hex1(1,1)=-Hex*m2(1,1);
Hex1(1,2)=-Hex*m2(1,2);
Hex1(1,3)=-Hex*m2(1,3);
Hex2(1,1)=-Hex*m1(1,1);
Hex2(1,2)=-Hex*m1(1,2);
Hex2(1,3)=-Hex*m1(1,3);
%anisotropy field
Hk1(1,1)=0;
Hk1(1,2)=0;
Hk1(1,3)=Hk*m1(1,3);   
Hk2(1,1)=0;
Hk2(1,2)=0;
Hk2(1,3)=Hk*m2(1,3);
%Heff effective field
    Heff1(1,1)=H0x+Hp(1,1)+Hex1(1,1)+Hk1(1,1);
    Heff1(1,2)=H0y+Hp(1,2)+Hex1(1,2)+Hk1(1,2);
    Heff1(1,3)=H0z+Hp(1,3)+Hex1(1,3)+Hk1(1,3);
    Heff2(1,1)=H0x+Hp(1,1)+Hex2(1,1)+Hk2(1,1);
    Heff2(1,2)=H0y+Hp(1,2)+Hex2(1,2)+Hk2(1,2);
    Heff2(1,3)=H0z+Hp(1,3)+Hex2(1,3)+Hk2(1,3);   
t(1,1)=0;
%LLG calculate the LLG equation
for i=1:1:n    
    t(i+1,1)=i*dt;

    rt1(i+1,1)=-1*gama*miu0*(m1(i,2)*Heff1(i,3)-m1(i,3)*Heff1(i,2))-miu0*alpha*gama*(m1(i,1)*m1(i,2)*Heff1(i,2)-m1(i,2)*m1(i,2)*Heff1(i,1)-m1(i,3)*m1(i,3)*Heff1(i,1)+m1(i,1)*m1(i,3)*Heff1(i,3))*at;
    rt1(i+1,2)=-1*gama*miu0*(m1(i,3)*Heff1(i,1)-m1(i,1)*Heff1(i,3))-miu0*alpha*gama*(m1(i,3)*m1(i,2)*Heff1(i,3)-m1(i,3)*m1(i,3)*Heff1(i,2)-m1(i,1)*m1(i,1)*Heff1(i,2)+m1(i,1)*m1(i,2)*Heff1(i,1))*at;
    rt1(i+1,3)=-1*gama*miu0*(m1(i,1)*Heff1(i,2)-m1(i,2)*Heff1(i,1))-miu0*alpha*gama*(m1(i,1)*m1(i,3)*Heff1(i,1)-m1(i,1)*m1(i,1)*Heff1(i,3)-m1(i,2)*m1(i,2)*Heff1(i,3)+m1(i,2)*m1(i,3)*Heff1(i,2))*at;
    m1(i+1,1)=m1(i,1)+dt*rt1(i+1,1);
    m1(i+1,2)=m1(i,2)+dt*rt1(i+1,2);
    m1(i+1,3)=m1(i,3)+dt*rt1(i+1,3);
    m1(i+1,4)=sqrt(m1(i+1,1)^2+m1(i+1,2)^2+m1(i+1,3)^2);
    m1(i+1,1)=m1(i+1,1)/m1(i+1,4);
    m1(i+1,2)=m1(i+1,2)/m1(i+1,4);
    m1(i+1,3)=m1(i+1,3)/m1(i+1,4);
    
    rt2(i+1,1)=-1*gama*miu0*(m2(i,2)*Heff2(i,3)-m2(i,3)*Heff2(i,2))-miu0*alpha*gama*(m2(i,1)*m2(i,2)*Heff2(i,2)-m2(i,2)*m2(i,2)*Heff2(i,1)-m2(i,3)*m2(i,3)*Heff2(i,1)+m2(i,1)*m2(i,3)*Heff2(i,3))*at;
    rt2(i+1,2)=-1*gama*miu0*(m2(i,3)*Heff2(i,1)-m2(i,1)*Heff2(i,3))-miu0*alpha*gama*(m2(i,3)*m2(i,2)*Heff2(i,3)-m2(i,3)*m2(i,3)*Heff2(i,2)-m2(i,1)*m2(i,1)*Heff2(i,2)+m2(i,1)*m2(i,2)*Heff2(i,1))*at;
    rt2(i+1,3)=-1*gama*miu0*(m2(i,1)*Heff2(i,2)-m2(i,2)*Heff2(i,1))-miu0*alpha*gama*(m2(i,1)*m2(i,3)*Heff2(i,1)-m2(i,1)*m2(i,1)*Heff2(i,3)-m2(i,2)*m2(i,2)*Heff2(i,3)+m2(i,2)*m2(i,3)*Heff2(i,2))*at;
    m2(i+1,1)=m2(i,1)+dt*rt2(i+1,1);
    m2(i+1,2)=m2(i,2)+dt*rt2(i+1,2);
    m2(i+1,3)=m2(i,3)+dt*rt2(i+1,3);
    m2(i+1,4)=sqrt(m2(i+1,1)^2+m2(i+1,2)^2+m2(i+1,3)^2);
    m2(i+1,1)=m2(i+1,1)/m2(i+1,4);
    m2(i+1,2)=m2(i+1,2)/m2(i+1,4);
    m2(i+1,3)=m2(i+1,3)/m2(i+1,4);
 %calculate the spin current   
    Is1(i+1,1)=m1(i+1,2)*rt1(i+1,3)-m1(i+1,3)*rt1(i+1,2);
    Is1(i+1,2)=m1(i+1,3)*rt1(i+1,1)-m1(i+1,1)*rt1(i+1,3);
    Is1(i+1,3)=m1(i+1,1)*rt1(i+1,2)-m1(i+1,2)*rt1(i+1,1);
    Is2(i+1,1)=m2(i+1,2)*rt2(i+1,3)-m2(i+1,3)*rt2(i+1,2);
    Is2(i+1,2)=m2(i+1,3)*rt2(i+1,1)-m2(i+1,1)*rt2(i+1,3);
    Is2(i+1,3)=m2(i+1,1)*rt2(i+1,2)-m2(i+1,2)*rt2(i+1,1);
    Is3(i+1,1)=m1(i+1,2)*rt2(i+1,3)-m1(i+1,3)*rt2(i+1,2);
    Is3(i+1,2)=m1(i+1,3)*rt2(i+1,1)-m1(i+1,1)*rt2(i+1,3);
    Is3(i+1,3)=m1(i+1,1)*rt2(i+1,2)-m1(i+1,2)*rt2(i+1,1);
    Is4(i+1,1)=m2(i+1,2)*rt1(i+1,3)-m2(i+1,3)*rt1(i+1,2);
    Is4(i+1,2)=m2(i+1,3)*rt1(i+1,1)-m2(i+1,1)*rt1(i+1,3);
    Is4(i+1,3)=m2(i+1,1)*rt1(i+1,2)-m2(i+1,2)*rt1(i+1,1);
    Is=(h/e)*(g11*Is1+g22*Is2+g12*Is3+g21*Is4);    
%Heff
    Hp(i+1,3)=0;
    Hp(i+1,2)=Hp0*sin(w*t(i+1,1)+thetay);
    Hp(i+1,1)=Hp0*cos(w*t(i+1,1)+thetax);
    Hex1(i+1,1)=-Hex*m2(i+1,1);
    Hex1(i+1,2)=-Hex*m2(i+1,2);
    Hex1(i+1,3)=-Hex*m2(i+1,3);
    Hk1(i+1,1)=0;
    Hk1(i+1,2)=0;
    Hk1(i+1,3)=Hk*m1(i+1,3);   
    Hex2(i+1,1)=-Hex*m1(i+1,1);
    Hex2(i+1,2)=-Hex*m1(i+1,2);
    Hex2(i+1,3)=-Hex*m1(i+1,3);
    Hk2(i+1,1)=0;
    Hk2(i+1,2)=0;
    Hk2(i+1,3)=Hk*m2(i+1,3);
    Heff1(i+1,1)=H0x+Hp(i+1,1)+Hex1(i+1,1)+Hk1(i+1,1);
    Heff1(i+1,2)=H0y+Hp(i+1,2)+Hex1(i+1,2)+Hk1(i+1,2);
    Heff1(i+1,3)=H0z+Hp(i+1,3)+Hex1(i+1,3)+Hk1(i+1,3);
    Heff2(i+1,1)=H0x+Hp(i+1,1)+Hex2(i+1,1)+Hk2(i+1,1);
    Heff2(i+1,2)=H0y+Hp(i+1,2)+Hex2(i+1,2)+Hk2(i+1,2);
    Heff2(i+1,3)=H0z+Hp(i+1,3)+Hex2(i+1,3)+Hk2(i+1,3);   
end
%extract the data every 10 points
tn=n/10;
for j=1:1:tn
    mt(j,1)=m1(j*10-9,1);
    mt(j,2)=m1(j*10-9,2);
    mt(j,3)=m1(j*10-9,3);
    mt(j,4)=m2(j*10-9,1);
    mt(j,5)=m2(j*10-9,2);
    mt(j,6)=m2(j*10-9,3);  
    Ist(j,1)=Is(j*10-9,1);
    Ist(j,2)=Is(j*10-9,2);
    Ist(j,3)=Is(j*10-9,3);  
    Is1t(j,1)=Is1(j*10-9,1);
    Is1t(j,2)=Is1(j*10-9,2);
    Is1t(j,3)=Is1(j*10-9,3); 
    Is2t(j,1)=Is2(j*10-9,1);
    Is2t(j,2)=Is2(j*10-9,2);
    Is2t(j,3)=Is2(j*10-9,3); 
    Is3t(j,1)=Is3(j*10-9,1);
    Is3t(j,2)=Is3(j*10-9,2);
    Is3t(j,3)=Is3(j*10-9,3); 
    Is4t(j,1)=Is4(j*10-9,1);
    Is4t(j,2)=Is4(j*10-9,2);
    Is4t(j,3)=Is4(j*10-9,3); 
end
T=1/fh;
Tn=T*1e-10/dt;
NI=floor(10*Tn);
Iselect(:,5)=Ist(tn-NI:tn,3);
Iselect(:,1)=(h/e)*Is1t(tn-NI:tn,3);
Iselect(:,2)=(h/e)*Is2t(tn-NI:tn,3);
Iselect(:,3)=(h/e)*Is3t(tn-NI:tn,3);
Iselect(:,4)=(h/e)*Is4t(tn-NI:tn,3);
% calculat the spin current along z axis
Isztot(thetaxn,1)=sum(Iselect(:,1))/10;
Isztot(thetaxn,2)=sum(Iselect(:,2))/10;
Isztot(thetaxn,3)=sum(Iselect(:,3))/10;
Isztot(thetaxn,4)=sum(Iselect(:,4))/10;
Isztot(thetaxn,5)=sum(Iselect(:,5))/10;
%save the data of magnetization and the spin current
strm1 = ['E:\AFM\m-h',int2str(hn),'thetax',int2str(thetaxn),'.xls'];
xlswrite (strm1,mt);
strm3 = ['E:\AFM\Is-h',int2str(hn),'thetax',int2str(thetaxn),'.xls'];
xlswrite (strm3,Iselect);
clear Iselect;
end
strm4 = ['E:\AFM\Istot-thetax-f',int2str(fh),'.xls'];
xlswrite (strm4,Isztot);
end